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We show that the parametrically driven nonlinear Schrodinger equation has wide classes of trav- 
elling soliton solutions, some of which are stable. For small driving strengths stable nonpropagating 
and moving solitons co-exist while strongly forced solitons can only be stable when moving suffi- 
ciently fast. 
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I. INTRODUCTION 

The parametrically driven damped nonlinear 
Schrodinger (NLS) equation, 



iipt + ipxx + 2|^)| ip — ip = hip* — i'jip, 



(1) 



was used to model the nonlinear Faraday resonance in 
a vertically oscillating water tank pl^j and the effect of 
phase-sensitive parametric amplifiers on solitons in opti- 
cal fibers §]. The same equation describes an easy-plane 
ferromagnet with a combination of a static and hf field 
in the easy plane ^,|| , and the planar weakly anisotropic 
XY model [j| . It also serves as a continuum limit for the 
parametrically driven Frenkel-Kontorova chain (an array 
of diffusively coupled pendula). The Frenkel-Kontorova 
system is regarded as a fairly realistic model of a num- 
ber of physical and biological systems and phenomena, 
including ladder networks of discrete Josephson junc- 
tions, charge-density-wave conductors, crystal disloca- 
tions in metals, DNA dynamics and proton conductivity 
in hydrogen-bonded chains . 

The undamped undriven nonlinear Schrodinger equa- 
tion exhibits soliton solutions which can travel with arbi- 
trary velocity and transport physical characteristics such 
as mass, momentum and energy. (For the sake of brevity, 
we are making use of the hydrodynamical interpretation 
of the equation here.) In equation ([!]), the second term in 
the right-hand side accounts for dissipative losses which 
occur in all physical systems. The dissipation has two 
visible effects on the soliton: it attenuates its speed and 
damps its amplitude. The main purpose of the intro- 
duction of the pumping (represented by the first term 
in the right-hand side) is to compensate for these losses. 
The parametric forcing is well known to be capable of 
counterbalancing the damping of the quiescent soliton's 
amplitude; a natural question now is whether it can sus- 
tain its motion with a nonzero velocity. 

In fact the existence of travelling solitons is a non- 
trivial matter even in the absence of damping. The 
driving term hip* in Eq.Q breaks the galilean invari- 
ance of the unperturbed nonlinear Schrodinger equa- 
tion and hence one cannot obtain a moving soliton sim- 
ply by boosting a static one. However the galilean or 
Lorentz symmetry is not a prerequisite for the existence 



of moving nonlinear waves. For example, dissipative sys- 
tems do not possess any symmetries of this kind but 
are well known to support stably propagating fronts and 
pulses whose velocities are determined by parameters of 
the model. In particular, travelling domain walls arise 
in the parametrically-driven Ginsburg-Landau equations 
(where the motion is due to nongradient terms) j8|. As 
far as solitons in Hamiltonian systems are concerned, the 
example of dark solitons in the nonlinear Schrodinger 
equations suggests that they have even greater mobility 
than dissipative fronts and pulses. Although in this case 
the galilean invariance is broken by the presence of the 
nonvanishing background, the dark solitons can propa- 
gate with arbitrary speeds bounded only by the velocity 
of sound waves p|-|ll"[] . 

A number of nonstationary regimes were reported in 
the water tank experiments, including the formation of 
oscillating soliton pairs but no steadily moving soli- 
tons were detected so far. On the other hand, numer- 
ical simulations of the undamped parametrically driven 
NLS equation [Eq .rtlJ) with 7 = 0] did exhibit travelling 
localised objects It has remained an open ques- 

tion whether these moving objects preserve their speed 
and amplitude, or attenuate and decay slowly due to the 
emission of the second-harmonic radiation. The aim of 
the present paper is to study the existence of steadily 
propagating solitons, and examine their stability. Here 
we are confining ourselves to the undamped situation rel- 
egating the analysis of the effect of damping to future 
publications. 

In addition to their role in transport phenomena, sta- 
bly moving solitons are also of interest as alternative 
attractors which may compete with (static or oscillat- 
ing) nonpropagating solutions. We will demonstrate that 
stable travelling solitons do exist in the parametrically 
driven nonlinear Schrodinger equation. Moreover, there 
are parameter ranges where moving solitons are stable 
whereas their quiescent counterparts are not. Unsta- 
ble solitons are not meaningless either; they arise as 
long-lived transients and intermediate states in spatio- 
temporal chaotic regimes. In this paper we will iden- 
tify oscillatory and translational instabilities of travelling 
solitons and simulate their nonlinear evolution near the 
transition curves. 

The structure of the paper is as follows. In section || 
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we derive an apriori bound for the existence domain of 
travelling solitons and introduce the linearised eigenvalue 
problem for their stability analysis. We also discuss some 
general properties of eigenvalues and eigenfunctions and 
formulate a simple criterion for the onset of the nonoscil- 
latory instability: dP/dV = 0, where V is the velocity 
of the steadily moving soliton, and P the associated mo- 
mentum. 



In section [II we present several explicit quiescent 
(V = 0) solutions, including a stationary complex of two 
solitons, and then derive the necessary condition for a 
static solution to be continuable to nonzero velocities. 
This condition is that the motionless solution should ei- 
ther not have any "free" parameters apart from the trans- 
lational shift, or, if there is an additional parameter z, the 
equation dP/dz = should be satisfied. Here P is the 
momentum of the motionless localised solution (which, 
contrary to one's mechanical intuition, is not necessarily 
equal to zero). There are three static solutions satisfying 
the above condition, two of which being the well-known 
constant-phase ip+ and ?/>_ solitons, respectively, while 
the third solution looks like a pulse with a bell-shaped 
modulus and twisted phase. 

The most important results of this work are contained 



in section IV where we report on the numerical contin- 



uation of various branches of solutions and their stabil- 
ity analysis. In agreement with the analytical predic- 
tions of the preceding section, we find that each of the 
above static solutions admits the continuation to nonzero 
V. The stability properties of travelling solitons result 
from an intricate interplay of two types of instability, the 
oscillatory and translational instability. In accordance 
with the conclusions of section ||, the numerical analy- 
sis of the linearised eigenvalues shows that the transition 
curves of the translational instability satisfy dP/dV = 0. 
One interesting conclusion of the stability analysis is 
that although quiescent solitons are unstable for driving 
strengths larger than h = 0.064, there are stable mov- 
ing solitons for any < h < 1. We discuss in detail 
the soliton's transformation as it is continued in V, pay- 
ing special attention to the dynamics of the associated 
linearised eigenvalues on the complex plane. Two dif- 
ferent scenarios of the transformation are identified, one 
occurring for small h and the other one for larger driving 
strength, and we also describe an interesting cross-over 
from one to another. 

Section ^ is devoted to the direct numerical simula- 
tions of the full time-dependent nonlinear Schrodingcr 
equation. We show that the evolution of both types of the 
soliton instability leads, as t — ► oo, to the same asymp- 



totic attractors. Finally, section VI contains concluding 
remarks and outlines some open problems. 



II. STEADILY TRAVELLING WAVES: 
EXISTENCE AND STABILITY 



A. Existence domain and integrals of motion 

We will confine ourselves to localised travelling waves 
of the simplest form, ip(X,t) = ip(X — Vt). Trans- 
forming to the moving frame, these correspond to time- 
independent soliton solutions of the equation 



iipt ~ iVip x + ip xx + 2\ip\ 2 ip -ip = hip* 



(2) 



where x — X — Vt. We will search for these static solu- 
tions by solving an ordinary differential equation 



iVip x + ip xx + 2|V>|V -ip = hip* 



(3) 



under the vanishing boundary conditions — > as 
|x| — » oo. Here h is always taken positive; negative h's 
can be recovered by the phase transformation ip — > iip. 

If ?/>(x) is a solution, then so is —ip(x). Next, it is 
straightforward to notice that if the function tp( x ) de- 
scribes a soliton travelling with the velocity V, the func- 
tion V'*( x ) describes a soliton moving with the velocity 
—V. We will try to restrict ourselves to positive V's 
wherever possible; we will only present negative velocities 
where this may help visualising how different branches of 
solutions are connected. Since the soliton moving with 
the velocity — V is given by tp(—x), the above observa- 
tion tells us that either tp*(x) = ±ip(—x) (that is, one of 
the real and imaginary part of the solution is even and 
the other one odd), or there are two solutions associated 
with the same V. (Here we are not making any difference 
between solutions which are different just in the overal 
sign.) In the latter case the solutions will not exhibit the 
ip*(x) = iip(-x) symmetry. 

Next, it is easy to show that solitons cannot travel 
faster than a certain limit speed. Indeed, as |x| — * oo, the 
soliton's asymptotic tail decays as ip(x) ~ e~ KX , where 



2k 2 = 2-V 2 ± 



(2 -V 2 ) 2 + 4(/i 2 -l) 



(4) 



Large driving strengths, h > 1, are of little interest to 
us as in this case the zero background, ip(x) = 0, is 
unstable with respect to continuous spectrum waves [pj. 
Therefore we are not going to discuss this case here. In 
the complementary region h < I, the complex struc- 
ture of k depends on the value of the velocity. When 
V 2 < 2 — 2\/l — h 2 , there are four real exponents; for 
2 - 2-s/l - h 2 < V 2 < c 2 , where 



= V2 + 2V1 -h 2 , 



(5) 



we have a quadruplet of complex k's, and finally, for 
V 2 > c 2 all four exponents are imaginary. Consequently, 
there can be no exponentially localised solitons travel- 
ling faster than c. Physically, c represents the minimum 
phase velocity of linear waves governed by Eq.(|l|), and 
our condition V < c is essentially an exclusion principle 
ruling out a resonance between solitons and linear waves. 

In the undamped case the parametrically driven NLS 
equation (0) conserves the momentum, 
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P= 2 / WZfl>-il>*1>*)dx, 



(6) 



and energy: 

E = ReJ (|^| 2 + |V| 2 - l^l 4 + dx. (7) 
Noting that stationary solutions satisfy 



2 L/.|2 i / 1 4 



+ /iReV' 



Eq.(^) can be rewritten as 

£ = 2Rc / (IV'I 2 - |^| 4 + /# 2 ) dx. (8) 



Containing no derivatives, this formula for energy of sta- 
tionary solutions has obvious advantages for the numer- 
ical implementation. 



B. Linearised eigenvalue problem 

In this paper we solve the equation ([}]) numerically and 
examine the stability of the resulting solutions by study- 
ing the associated eigenvalue problem. This eigenvalue 
problem arises by assuming a small perturbation of the 
form 

5ip(x, t) = y(x)e xt , y(x) = Su(x) + iSv(x). 
Substituting into (^) gives 

HY = XJY, (9) 
where the hermitean operator TL has the form 



H = I{-d 2 x + 1) + VJd x + 



h-6u 2 - 2i 
—4uv 



—4uv 
—h — 6v 2 — 2u 2 



the matrix J is given by 



J = 



(10) 



(11) 



and the column-vector Y(x) — (Rey,Imy) T = (Su,6v) T . 
In Eq.(^) / is the identity matrix, and we have decom- 
posed the stationary solution as ip(x) = u{x) + iv(x). 

For symmetric solutions satisfying ip*(x) — ±ip(—x) 
eigenvalues will always come in (A, — A)-pairs. This 
follows from the fact that for these solutions chang- 
ing x — > —a; in the operator ( |l0| ) amounts to chang- 
ing the sign of its off-diagonal elements, and hence if 
(5u(x), Sv(x)) T is an eigenfunction associated with an 
eigenvalue A, the column (Su(— x), — Sv(—x)) T will serve 
as an eigenfunction associated with an eigenvalue —A. 
As far as the zero eigenvalue is concerned, it will have 



a twin with the eigenfunction (Su(— x), — 6v(— x)) T un- 
less its eigenfunction y = 6u + iSv satisfies the symmetry 
y*(x) = e Iip y(—x), where if = const. 

To complete the discussion of the spectrum structure, 
we mention that there are two branches of the continuous 
spectrum lying on the imaginary axis of A: A = iuii^ik), 
where 

wi, 2 (fc) = Vk± V(fc 2 + l) 2 - h 2 , 

and — oo < k < oo. (We are still assuming h < 1). In 
the region V 2 < c 2 which is of interest to us, the con- 
tinuous spectrum has a gap: u>i(k) > ujq, u>2(k) < —ojq, 
where loq > 0. This gap can harbour discrete eigenvalues 
representing stable oscillation modes. 



C. Translational (nonoscillatory) instabilities 

The aim of this subsection is to demonstrate that a 
pair of pure imaginary eigenvalues can collide at A = 
and move onto the real axis only at the point where 
dP/dV — 0. This criterion is known in the context of 
dark solitons of the undriven nonlinear Schrodinger equa- 
tions; see pel , [ fyj . Here we simply adapt the proof given 
in pp[ | to the case of the equation with the parametric 
forcing. An important assumption that we make here, 
is that the solution whose stability is being examined, 
does not have any free parameters apart from the trivial 
translation parameter, xo- 

First of all we need to make a remark on the integrable 
case, h = 0. In this case solutions of the ODE (S) can be 
obtained from a quiescent soliton of Eq.(|l|) by a Galilei 
transformation: 



fj,(x) = e l{v/2)x AsechAx, 



(12) 



where A = - V 2 /^. For h = and any \V\ < 2, 
the linearised operator H. has four zero eigenvalues asso- 
ciated with two eigenvectors. One of these eigenvectors 
originates from the translation symmetry and the other 
one results from the phase invariance of Eq. (|^) . The term 
hip* breaks the phase invariance and hence as h is in- 
creased from zero, one pair of eigenvalues (A, —A) moves 
away from the origin on the complex plane. As h and 
V are further varied, a pair of eigenvalues may return 
to the origin. If the solution of Eq.([5|) at the point of 
their return is a member of a family parametrised by two 
free parameters, we will have, again, four zero eigenvalues 
with two eigenfunctions. (The eigenfunctions are simply 
derivatives of the solution with respect to the free param- 
eters.) Our analysis will not be applicable in this case, 
and the equality dP/dV = does not have to be valid 
at the return point. (We will come across this type of a 
situation in section IV D| below.) However, a more com- 
mon situation is when the solution at the return point is 
a member of a one-parameter family. We will show that 
in this case the relation dP/dV = does have to be in 
place. 
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Let us denote V c the velocity for which the eigenvalue 
of the operator @-© vanishes. We can develop the 
solution i/)(V; x) in powers of e = V — V c : 

ip(V; x) = f {x) + ei/>t(x) + e 2 ^ 2 {x) + 

where tpo = if)(V c ]x). Accordingly, the operator Tl ex- 
pands as H = Ho + eHi + e 2 H 2 + .... If the eigenvalue 
A moves from imaginary to the real axis, it is natural to 
assume that it admits an expansion of the form 

A = e^Ai + e 3 / 2 A 3 + e 5 / 2 A 5 + . . . . (13) 

The associated eigenfunction is then developed as 

Y{x) = Y (x) + e 1/2 Y x {x) + eY 2 {x) + .... (14) 

When e = 0, we have HqYo = 0, i.e. Yq is a null eigen- 
vector at the bifurcation point V = V c . Since we have 
assumed that tp{V c , x) is a member of a one-parameter 
family of solutions, the operator Ho has only one null 
eigenvector, and we have to identify Y"o = ^J (x). Here 
is a column- vector formed by the real and imaginary 
part of the soliton ip : ^ = (u ,vq) t . The prime indi- 
cates differentiation with respect to x. 

Next, setting the coefficient of e 1 / 2 to zero yields 

H.a Yi = Ai JYq. 

Comparing this to the equation 

<9* 



Ho 9V 



v=v c 



which arises from the differentiation of Eq.(||) with re- 
spect to V, we get 



M,) = -A 1 - 



v=v c 

T 



(In the above equations "J" = (it, v) .) The coefficient of 
e 1 produces 

H Y 2 =X 1 JY 1 -H 1 Y , 

which has bounded solutions if the right-hand side is or- 
thogonal to the null eigenvector of Ho'- 



Ai / YoJYtdx - / Y HiY dx = 



(15) 



The second term in equation (|15|) is readily shown to van- 
ish — one only needs to expand the identity 'WS' = 
in e. (The coefficient of e 1 gives Hi^' = -Ho^[. 
Taking the scalar product with yields the required 
J ^'Mi^'^dx — 0.) On the other hand, the first term in 
Eq.@ is equal to {\\/2)dP/dV. Consequently, Eq.@ 
gives either dP/dV — or Ai = 0. If we assume that 
Ai = 0, we will not be able to conclude that dP/dV = 
at this order of the expansion. However, the order e 2 



will then give us X^dP/dV — 0, which implies either 
dP/dV — or X 2 = 0. Proceeding by a similar token we 
will eventually arrive at the equation dP/ dV = at some 
order e™ where n is such that A„ ^ 0. (Alternatively, we 
will have to conclude that all A„ = and hence we are 
dealing with a symmetry eigenvalue which is equal to 
zero for all V.) 

Thus a pair of real or pure imaginary eigenvalues of 
the same magnitude and opposite sign, can only collide 
for the value of V which satisfies dP/dV — 0. Here we 
wish to re-emphasise that we have obtained this conclu- 
sion under the assumption that the geometric multiplic- 
ity of the zero eigenvalue is not increased at the point 
of collision. A simple example when this assumption is 
not valid, is furnished by the case h = 0. In this case 
the momentum corresponding to the soliton (|l^ ) is given 
by P = V 'y/l — V 2 /A. Although P has a maximum for 
V = s/2, the stability properties of the undriven soliton 
do not change at this point. The reason is that for each V 
the operator H has two null eigenvectors in this case, and 
hence we cannot make the identification Y = ^' . (In- 
stead, Yq will be a linear combination of two zero modes.) 
Consequently, the above proof becomes invalid in this 
case. 

Finally, one can easily check that the above result does 
not really depend on how the eigenvalue A expands in 
powers of e. We assumed that the expansion (|l^) starts 
with terms of order e 1 / 2 . This assumption is natural and 
supported by the numerical evidence; however, even if we 
had postulated the expansion starting with terms of order 
e 1 / 4 , e 1 / 3 or say, e, we would have still arrived at the same 
necessary condition for the zero crossing: dP/dV = 0. 



III. QUIESCENT SOLUTIONS (V = 0) 

A. The "twist" soliton 

In order to continue in V it is useful to have some start- 
ing solutions for V = 0. Two such stationary solutions 
are given by 



ip+(x) = A + sech(A + x), 
ip-(x) — iAsech(A-x), 



(16a) 
(16b) 



where A\ = 1 ± h. The soliton ■0- is unstable with re- 
spect to a nonoscillatory mode for all ft, ||. The tp + is 
stable for h < ho — 0.063596 but developes an oscillatory 
instability as h is increased beyond ho [p|p^] . 

In this section we will produce several more explicit 
solutions of the undamped, parametrically driven NLS 
equation (Q). Writing ip = u + iv, the stationary equa- 
tion (|^) transforms into the system 



u X x — u — hu + 2u(u + v ) = 0, 
Vxx — v + hv + 2v(u 2 + v 2 ) = 0. 



(17) 
(18) 
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We can try to find explicit solutions of this system by 
imposing some plausible reductions, for instance by iden- 
tifying u 2 + v 2 with a function of u and its derivatives: 



u x , Uxx and so on. In this case Eq.(|17]) is an equation 
for u only, while Eq.([l8|) becomes a linear equation with 
variable coefficients. The simplest choice u 2 + v 2 = Cu 2 
leads to the ip+ and ip- solitons Another simple 

possibility is to require that 



Cu, C — const; 



(19) 



this converts the first equation into the stationary KdV 
with the well-known localised solution 



3 1 



c 



— sech 2 £, 



The second equation has now the form of an eigenvalue 
problem for the Poschl- Teller potential: 



-91 



1 - 6sech 2 £)« = Ev, 



The operator 

2^ 



(20) 
has 



where E = 1 + A(h - l)/(h + 1 

two localized eigenf unctions, vq = asecli^£ associated 
with an eigenvalue Eq = —3, and v\ = asech£tanh£ 
with an eigenvalue E\ = 0. The first cigcnfunction can 
satisfy the constraint ( |19|) for no a while the second one 
will satisfy it if we set C 2 — a 2 = |(1 + h). 



itisfy it if we set C 2 = a 2 = |(1 + h). 
Noticing that E = corresponds to h = |, 



we con- 



clude that for h 
form 



ut 



^ sech 2 ^; 



there is an explicit solution of the 



(21) 



vt = ±\/ - sech£ tanh^, 



where £ = y ^x. Similarly to the solitons ip+ an d ip—, 
the modulus of the above solution is bell-shaped, but, 
unlike the constant phase of the ip+ and ijj- solitons, the 
phase of this solution varies. (In the case of the positive 
sign in ( |2l| ) the phase grows from —ir/2 at x — — oo to 
7r/2 at x = oo.) The solution looks like a pulse twisted 
by 180° in the (u, w)-plane. For this reason we will be 
referring to Eq.(|2l|) as the "twist" soliton. 
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= 2(Si+S 2 -2/3) 



D = 1 + e"" L + er" + e~ 
dx = A_{x-z), 6 2 = A + {x + z); 



B. The twist soliton as a bound state 

The system (|l7|)-(|l8|) appeared previously as a station- 
ary system governing light pulses in a birefringent opti- 
cal fiber. Using Hirota's approach, Tratnik and Sipe fl3] | 
have obtained the following exact solution to eqs.p),(|]j 



tp — ip(z; x) = u + iv, 



o 2(0i-/3) 



u = 2A+e 02 D- 1 (1 + e 
v = 2A^D- 1 (l-e 2 ^-V 



(22a) 
(22b) 

(22c) 



where 



the constant 



A_ 



Aa 



> 0, 



(23) 



the amplitudes A± are as in ©: A± = VT±h, and z 
is a real parameter which can take arbitrary values. The 
solution ( p2| ) with z = 10 and —10 is plotted in Fig.l, (a) 
and (c). As is clear from the figure, for sufficiently large 
\z\ the solution represents a complex of two solitons, 
and ip-, with the separation approximately equal to 2\z\. 
It is perhaps worth emphasizing here that the parameter 
z has nothing to do with shifting the solution as a whole, 
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x — > x — xq, which is possible due to the translation in- 
variance of Eqs.(0)-@. (This overall shift parameter, xo, 
is disregarded in equations (21 )-(p2|) and in the remain- 
der of this text). The parameter z is nontrivial in the 
sense that the shape of the solution depends on z. 

It is straightforward to verify that our "twist" solu- 
tion is a particular case of Eq.(p2|) with h = |. Indeed, 

choosing z = — g-y^f In 3 in Eqs.(|22|), we obtain the soli- 
ton (the one with the negative sign), centered at 
x = — 3z. Since the twist is a symmetric solution with 
secft-shaped modulus, it would be difficult to interpret it 
as a bound state of the ip+ and tp— without embedding 
it into a broader family of solitonic complexes. In fact, a 
similar symmetric solution exists for any h, not only for 
h = |. To see this, we notice a simple relation between 
two solutions of the form ( |22] ) — one with the parameter 
value z = C + £ and the other one with z = ( — £: 

i>(t + br]-y) = r((-Z;v + y)- (24) 

Here f and r\ are defined by the driving strength h: 
(3(1 1 



c 



2\A_ Aj 



< 



and 



V 



a (± 

2 \A- 



A A 



(25) 



(26) 



while £ and y can take arbitrary values. The relation ( |24| ) 
implies that the solution (|2^) with z — £ is symmetric 
about the point x = rj: 



^(Cv-y) =ip*(C;'n + y)- 



(27) 



That is, the real part of this solution is even and imagi- 
nary part odd with respect to x = ?y: 

u(v-y) =u(v + y), v(r]-y) = -v(r) + y). 

(See Fig.l (b).) This particular representative of the fam- 
ily (|2^) will play a special role in what follows. For the 
ease of reference we are retaining the name "twist" for 
this symmetric solution — for all h. 



C. The moving soliton bifurcation 

Suppose the equation (^) has a one-parameter family 
of quiescent solutions ij)(z\x). Here z can be any non- 
trivial parameter; the only requirement is that z should 
not be mst an overall shift in x. One such family is given 
by Eq.(2^) and there can also be other families for which 
-0 is not available explicitly. We will show in this sec- 
tion that in order for a solution with some z = z to be 
continuable to nonzero V, the corresponding momentum 
integral should satisfy 



OP 

dz 



0. 



(28) 



Let us assume that Eq.(||) with V^O has a solution 
ip(x), and that this solution is an analytic function of V 
in some neighbourhood of V = 0. Then we can expand 
it in the Taylor series 

xl>(x) = Mx) + Vih(x) + V 2 ip 2 (x) + (29) 

where ipo(x) = tp(zo;x) = uq + iv^ is some representa- 
tive of the family of "motionless" solutions ip(z;x) with 
the parameter value zq. Substituting (|2^) into (g) and 
equating coefficients of like powers of V, we get, at the 
order V . 



n 



Ul 



= Jd x 



Uq 



(30) 



Here u\ + iv\ — ipi and the operator 7i is given by 
Eq. (|l0|) . The equation ( p0| ) is solvable in the class of 
square integrable functions if the vector in the right-hand 
side is orthogonal to all homogeneous solutions, i.e. to all 
null eigenvectors of the operator Ti. Since there is a fam- 
ily of "motionless" solutions parametrized by z and by 
an arbitrary spatial shift xq (which we have disregarded 
so far), the operator Ti. has two zero modes. One is the 
translation mode d x ipo — d x (uo+ivo)', the corresponding 
solvability condition is trivially satisfied: 



d x (u ,v )Jd x 



u 



dx = 0. 



The other zero mode is given by the derivative d z ipo = 
d z uo + id z VQ. The associated solvability condition reads 

IdP 
'2~dz~ 7 



= J {d z u ,d z v )Jd x 



Uq 



dx 



where P is the momentum integral (g). Consequently, a 
solution with nonzero V can only detach from the V = 
branch at the point where dP/dz = 0. 

Coming back to our explicit solutions, the ip+ and tp- 
solitons do not have any free parameters apart from the 
trivial position shift. Consequently, both solutions are 
continuable to nonzero V. Next, we have a family of 
solitonic complexes (|22|) with a nontrivial parameter z. 
As one can easily check, the momentum of the complex 
( p2| ) as a function of z has a single minimum for some 
finite z = zq and tends to zero as z — > ±oo. To find Zq, 
we notice that the relation ( |2^ ) implies 

p(C + = P(C - 0- 

This means that the function P(z) is even with respect 
to the point z = ( and therefore, ( is the point of the 
minimum: zq = £. Thus, the only representative of the 
family of the two-soliton complexes (22) that can be con- 
tinued to nonzero V, is our twist soliton, ip(C, x). (To be 
more precise, there are two twist solutions, one with pos- 
itive and the other one with negative momentum. This 
is related to the fact that when V = 0, we can generate 
new solutions to the system (p7|)-(p^) by changing the 
sign of just one component u or v.) 



6 



IV. BIFURCATION DIAGRAM 



A. The travelling ip- soliton 



We used a predictor-corrector continuation algorithm 
with a fourth-order Newtonian solver to continue solu- 



tions of equation 
mentum integral (t 



in V. Since derivatives of the mo- 
J) determine stability and branching 
properties of solutions, the momentum was our natural 
choice for the bifurcation measure. Eq.(||) was solved 
under the vanishing boundary conditions tp (±L/2) = 0. 
We used L = 200 (except the cases where we had to ex- 
tend the interval to account for slow decay of solutions) 
and the discretisation stepsize Ax = 0.005. The eigen- 
value problem (H) was solved on the interval (—50,50). 
Here we utilised the Fourier method, typically with 600 
harmonics. 



We start our description with the branch departing 
from the quiescent soliton t/>_. For every h this branch 
continues all way to V = c, where c is the minimum phase 
velocity of linear waves given by Eq.(|). (See Fig.2). As 
V — > c, the decay rate of ip(x) decreases and the soliton 
merges with the zero solution, with the momentum P 
tending to zero. We plotted P(V) for h = 0.1, 0.3, 0.5, 0.7 
and 0.9 in Fig.2. For technical reasons we could not con- 
nect the curve P(V) to zero although we were able to 
approach the value V = c as close as the fourth digit 
after the decimal point. (The problem is that since the 
decay rate of the solution decreases, one has to increase 
the length of the integration interval — and this cannot 
be done indefinitely.) The only curve which is connected 
to zero in Fig.2, is the one for the undriven case h = 0. 
In this case we enjoy an explicit solution ( O ) with the 
momentum P = V 'yl — V 2 /4. 



1 .8 



^,h=0.9 



*K-,h=0.2802 



W + ,h=0.27 



*K-,h=0.1 

^,h=0.05 
h=0 




Fig.2 The momentum of the ip+ and ip- solitons as a function of their velocities. Thick respectively thin lines depict 
stable respectively unstable solutions. Note that as h — > 0, the stability domain of the ?/>- tends to y2 < V < 2 and that 
of the ip + to < V < V2. The whole of the h = branch is stable. 



For each h the momentum of the soliton has a single 
maximum on this branch, at V — V c (Fig.2). To the left 
of V c the linearised operator (|9|)-(^0|) has a pair of real 
eigenvalues ±A and consequently, the soliton which 
is well known to be unstable for V = ||, remains un- 
stable for small nonzero velocities. As V approaches V c , 



the two eigenvalues converge at the origin on the com- 
plex plane, with the associated eigenfunctions tending to 
the translation mode ^>' (x). Increasing V past V c , the 
eigenvalues move onto the imaginary axis and hence the 
-0- soliton becomes stable for sufficiently large velocities 
(where dP/dV < 0). This change of stability proper- 
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ties is i n exa ct agreement with the scenario described in 
section 



IIC 



B. The travelling lp+ soliton; h < 0.25 



In agreement with conclusions of section [II C , we have 
found that the soliton ?/> + is also continuable to nonzero 
velocities. Unlike the -^--branch where the final product 
of the continuation process does not depend on the value 
of h, the transformation scenario of the soliton ip+ i s 
different for different h. For h < 0.28 the fate of the 
soliton ip + is similar to that of the ip-. As V — > c, 
the soliton develops oscillations on its tails; the width 
of the resulting oscillatory "wavepacket" grows and the 
amplitude decreases, until the solution becomes equal to 
zero everywhere. The momentum P(V) tends to zero as 
V — > c and has a single maximum at some V — V c . Sta- 
bility properties of the tp+ soliton depend on whether h 
is smaller than 0.064, lies between 0.064 and 0.25, or is 
greater than 0.25. 

Let, first, h < 0.25. As V is increased to V c , a pair 
of imaginary eigenvalues ±A collides at the origin on the 
complex plane and moves onto the real axis. This is in 
contrast to the soliton, where two real eigenvalues 
converged at the origin as V was increased. For V > V c 
(i.e. where the slope dP/dV is negative), the ip + soli- 
ton becomes unstable. (Note that the -0- was unstable 
for positive slopes dP/dV.) As in the ip- case, this is 
the translational, or nonoscillatory, instability. Since, as 
we checked, the eigenfunctions associated with the col- 
liding imaginary eigenvalues tend to ^f' (x) as V — > V c , 



this instability is of the type analysed in section II C 



The numerically detected instability boundary, defined 
by the condition dP/dV = 0, is in exact correspondence 
with our analytical predictions. 

In addition, for some V the ip+ soliton exhibits the 
oscillatory instability. (The oscillatory instability sets in 
when four eigenvalues collide, pairwise, on the imaginary 
axis and emerge into the complex plane.) For h < 0.064 
the oscillatory instability does not arise for any V and 
hence the whole range < V < V c is stable. (See the 
h = 0.05 curve in Fig.2). For 0.064 < h < 0.25 the 
oscillatory instability occurs for V zero and small, but 
as V is increased, both imaginary and real parts of the 
"unstable" eigenvalues decay, with the real parts decay- 
ing faster. Eventually the eigenvalues ±A, ±A* converge, 
pairwise, on the imaginary axis and the soliton stabilizes. 
Increasing V still further, two of the resulting imaginary 
eigenvalues, A and —A, start approaching each other. At 
V = V c where dP/dV = 0, they collide and move onto 
the real axis. The soliton looses its stability once again 
— this time to a nonoscillatory mode. This scenario is 
exemplified by the curves h = 0.1 and h = 0.15 in Fig.2. 

An interesting phenomenon occurs in the undriven 
case, h = 0. As we mentioned, in this case the sta- 
tionary solution of Eq.(0) is available in explicit form, 



Eq.(|12(). We also noted that the associated momentum 
P = VW1 — V 2 /4 has a maximum for finite V = 
On the other hand, in any small-ft. neighbourhood of this 
solution there are solitons "0+ and ip-, of which tp + is 
stable only for dP/dV > (i.e. for V < V2) and tp- is 
stable only for dP/dV < (i.e. for V > V2). How can 
this be reconciled with the fact that the solution ( |l2| ) is 
stable for all V for which it exists (V < 2)1 

The answer is related to the behaviour of the zero 
eigenvalue associated with the phase invariance of Eq. (m) 
with h — 0. As h increases from zero, this eigenvalue 
moves away from the origin. On the tp+ branch, it passes 
onto the imaginary axis, while on the ip- branch, the 
eigenvalue moves onto the real axis instead. 



C. The travelling i/>+ soliton; h > 0.25 

Let now 0.25 < h < 0.28, and assume we are moving 
along the ip+ branch in the direction of larger V. For 
small V we have a quadruplet of complex eigenvalues 
±A, ±A* implying the oscillatory instability. As V is in- 
creased, both imaginary and real parts decay — as in the 
h < 0.25 case. However, this time imaginary parts decay 
faster than the real parts, and the two pairs of eigen- 
values converge on the real axis. For velocities above 
this point the oscillatory instability is replaced by the 
translational instability. As V is increased further, one 
pair of the newly born real eigenvalues grows in absolute 
value whereas the other pair decreases in magnitude. At 
the point V — V c where P(V) reaches its maximum, the 
latter pair converges at the origin and moves onto the 
imaginary axis. (This does not render the soliton stable 
though, as the other pair remains on the real axis.) This 
scenario is exemplified by the curve h = 0.27 in Fig.2. 

Next, let h be greater than 0.28. For these h the branch 
P(V) emanating from the origin, turns back at some 
V — Vmax (Fig.2), with the derivative dP/dV remaining 
strictly positive for all V < V max . Below we will describe 
the transformation this solution undergoes when contin- 
ued beyond the "turning point" , while here we only wish 
to emphasize that no new zero eigenvalues can appear at 
this point. The reason is that V = V max is a bifurcation 
point of solutions of the ordinary differential equation (|^) 
but not of the partial differential equation (Q)-(||). (In 
other words, V is an "internal" parameter characterising 
the solution and not an "external" control parameter.) 
Indeed, the soliton is a member of a two-parameter (xq 
and V) family of solutions of Eq. ([[])- (||) and hence for 
any V there are two zero eigenvalues in the spectrum of 
the linearised operator Ti. Consequently, although being 
a turning point for the ODE (||), the value V = V max 
is in no way special as far as the PDE (0)-(||) and its 
linearisation are concerned. No changes of the soliton's 
stability properties occur at this velocity. 

How does one type of behaviour of the curve P(V), 
occurring for h > 0.28, replace the other one, arising 
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for h < 0.28? We scanned the interval 0.28000 < h < 
0.28020 and discovered a tiny region of transitional be- 
havior, around h — 0.28005. For this h, P(V) grows 
until it reaches a maximum at V c = 1.051 and then 
starts decreasing, as in the case of h < 0.28. However, 
the curve does not decay all way to P = as would 
be the case for h < 0.28, but reaches a minimum at 
V cc — 1.0563. After that, the momentum starts growing, 
and, at V max — 1.0565, the curve P(V) turns back - 
as for h > 0.28! To give an idea of how small this win- 
dow of transitional behaviour is, it suffices to say that for 
h = 0.28000 the momentum P(V) tends to as V -> c, 
whereas for h as close as 0.28010, the curve P(V) has a 
"turning point" . 




0.2 0.4 0.6 0.8 1 

Fig. 3 Stability diagram for the ip + and solitons on 
the (h, V)-plane. In the region marked "stable" one of the 
two one-soliton solutions is stable whereas the other one is 
not. Across the solid line, the corresponding soliton looses 
its stability to an oscillatory or monotonically growing mode. 
No solitons exist in the dashed region. 

What happens to the ip + soliton with h > 0.28 (more 
precisely, with h > 0.28010) as we continue it beyond 
the turning point? Fig. 4 (a) shows the momentum as a 
function of V. The point of intersection with the vertical 
axis V = corresponds to the twist solution (Eq.(|2^) 
with z = C given by Eq.(p5|).) On the diagram Fig. 4, it 
is marked as ipT- As we mentioned before, the (V = 0)- 
twist is a representative of a two-parameter family of sta- 
tionary solutions of Eq.(|^). Consequently, there should 
be four zero eigenvalues in the spectrum of the operator 
TL in this case, with two linearly independent eigenfunc- 
tions given by d x ^(z]x)\ z=( . and d z ^{z; x)\ z _^. (Here 
^f(z;x) is a two-component vector formed by the real 
and imaginary parts of Eq.(^).) Numerically, we ob- 
served that as we approach the (V — 0)-twist from the 
direction of positive V, a pair of opposite eigenvalues 
converges at the origin on the complex plane. The curve 
P(V) does not have an extremum at this point and this 
may seem to be in contradiction with predictions of sec- 
tion |[I C| . The paradox is resolved as soon as one recalls 
that the extremality condition dP/dV — was derived 
under the assumption that there is only one eigenvec- 
tor associated with the zero eigenvalue whereas we have 
two linearly independent null eigenvectors in the case at 



hand. 

As we continue further into the region V < 0, the twist 
gives rise to a variety of multisoliton complexes; we shall 
describe them in the next subsection. Here we will re- 
strict ourselves to the region V > where this branch can 
still be regarded as a branch of one-soliton solutions. Al- 
though these solutions undergo similar transformations 
for all h in the interval (0.28, 1), there are a few differ- 
ences with regard to the trajectories of eigenvalues on the 
complex plane. One difference worth to be mentioned 
is that for the driving strength h = 0.3 and larger h, 
the quadruplet of complex A persists on the entire upper 
branch of P(V) (i.e. for all V > 0). This is in contrast to 
the case of 0.25 < h < 0.28, where the complex quadru- 
plet converges on the real axis. Near the left end of the 
interval 0.28 < h < 1 (e.g. for h = 0.2802), we have 
an intermediate pattern. Similarly to the case h < 0.28, 
here the complex quadruplet converges on the real axis 
somewhere on the lower branch of P(V) (i.e. before the 
turning point), but as we move onto the upper branch, 
the two emerging real pairs reunite quickly and the com- 
plex quadruplet reappears. 

Next, as we know, there are only two points where a 
pair of eigenvalues can pass from the real onto the imag- 
inary axis, or vice versa. One point is V = V c where 
dP/dV = 0, and the other one is V — 0. Therefore 
the dynamics of eigenvalues depends on which of the two 
points comes first, or, equivalently, whether the upper 
branch of P(V) has the maximum for positive or nega- 
tive V. For smaller values of h in the interval (0.28, 1) 
(e.g. h = 0.3), where V c > 0, two imaginary eigenval- 
ues move to the real axis at V = V c . These imaginary 
eigenvalues have detached from the continuous spectrum 
somewhere before the turning point (i.e. on the lower 
branch of P(V).) The two newly born real eigenval- 
ues first diverge from the origin but then reverse and, 
at V — 0, move back onto the imaginary axis. For larger 
h (e.g. h = 0.7), where V c < 0, the pattern is different. 
Here, the two imaginary eigenvalues become real not at 
the point V c but at V = 0. Subsequently, as we con- 
tinue the branch to negative velocities, another pair of 
imaginary eigenvalues detaches from the continuum and 
at the point V c < two (imaginary or real) eigenvalues 
pass through the origin. 

Fig. 3 shows the stability diagram of the ip+ and ip- 
solitons on the (h, V)-p\a,ne. For the ip + soliton, the 
range of stable velocities approaches < V < \/2 
as h — > 0, while the stability range of ip- tends to 
V2 < V < 2. Finally, the domain of stability in the h = 
case is the union of the above two ranges: < V < 2. 



D. Other branches; h > 0.28 

As we continue it to negative velocities, the twist (we 
are using this name here for V ^ 0-dcformations of the 
quiescent twist solution) gradually transforms into a com- 







plex of two twists (plotted in Fig. 5(a)). A further con- 
tinuation of this branch takes us, via several "turning 
points" , to a solution that can be interpreted as an as- 
sociation of the twist and two ip- solitons of opposite 
polarities (denoted ip(-T-))- This solution is depicted in 
Fig.5 (b). 
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Fig.l (a): Tr)e full bifurcation diagram for h p= Ofl. All 
branches shown in this figure are unstable. (For h = 0.7 the 
soliton is stable at large velocities but the stable region 
lies beyond the frame of this figure — see Fig. 2.) More so- 
lution branches can be obtained by the reflection V — > —V, 
P — ► —P. (b): The corresponding E{P) dependence. The 
stable branch of i/j- is depicted by a thick line at the bottom 
of the figure. 

Fig. 4 (b) shows the energy of different branches as cal- 
culated by Eq.(||). We have eliminated the dependence 
on the soliton's velocity between P(V) and E(V) to ob- 
tain E as a function of P. The purpose of this "Leg- 
endre transformation" is the following. The stationary 
equation (|^) can be regarded as a condition that the en- 
ergy (@) be stationary under the fixed momentum (0): 
(8E)p = or, equivalently, S(E — VP) = 0, where V is 
the Lagrange multiplier. Consequently, a steadily travel- 
ling soliton satisfies the relation of the Hamilton mechan- 
ics, V = dE/dP, and its velocity is given by the slope of 
the corresponding E(P) branch in Fig. 4 (b). The relation 
SE = V8P implies that the curves E(V) and P(V) have 
extrema at the same set of points V = V C which will ap- 
pear as cusps on the E{P) plot. Since the function P(V) 
has extrema precisely at points where linearised eigen- 
values move from the imaginary to real axis, branches of 
the E(P) curve separated by the cusps may have different 



stability properties. One such change of stabilities does 
indeed occur in Fig. 4 (b) where the i/j_-curve is seen to 
be partitioned into stable and unstable branch. Solitons 
on the stable branch have lower energies than unstable 
solitons with the same values of the momentum. 
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Fig.5 (a): The V'(tt) complex, (b): The ip(-T-) solu- 
tion, (c): A complex of two twist solitons arising from the 
continuation of the "0(++) bound state. (Note the differ- 
ence from the other two-twist complex shown in (a).) Solid 
line: real part; dashed line: imaginary part 



Another branch emanating from the origin in Fig. 4 (a), 
is a bound state of two solitons ip+ . This solution was not 
obtained by the continuation from V — as Fig. 4 may 
seem to be suggesting. Instead, we fixed a nonzero V and 
continued in h from the value h = 0.05 where the com- 
plex ^(++) arises f rom t he ^-continuation of the twist 
soliton (see section IV E). Omitting details of this pro- 
cedure, we start the description of the resulting branch 
at some point (V, P) away from the origin. As we ap- 
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proach the origin from this point, the separation between 
the solitons ip+ in the complex t/'(++) rapidly increases 
so that the field values between the two solitons become 
exponentially small. For example, for h = 0.7, the (nu- 
merically calculated) separation at the point V = was 
equal to z ps 21. The value of at the point on the 
x-axis, equally distanced from the left and right soliton, 
was of order 10 -6 . Consequently the nonlinear term in 
the equation (||) becomes negligible away from the soli- 
tons' core and, in spite of an extremely small value of the 
residual that we used in our numerical algorithm (10~ 10 ), 
we were unable to distinguish between a genuine bound 
state and a linear superposition of two distant solitons. 
We conjecture that the complex ipu-+) exists all way to 
V = but as V — » 0, the intersoliton separation z — > oo. 
Another indication to this effect is that as V — > 0, the 
imaginary part of the solution tends to zero, rapidly and 
uniformly. For example, the imaginary part of the above- 
mentioned numerical solution with the real part between 
the solitons > 10~ 6 , was smaller than 10~ 13 for all x. 
Since the only pure real solution that exists for V = is 
the (single) soliton ip+, the V — > limit of the 4>(++) com- 
plex should be an infinitely separated pair of the ip+ , s. 

If we, conversely, continue our solution away from the 
origin, the curve P(V) turns left at some V and the com- 
plex tp(++) transforms into what can be interpreted as a 
bound state of two twists (denoted V'(tt) m Fig. 4.) This 
solution is depicted in Fig. 5(c). As V — > 0, the momen- 
tum of this bound state tends to zero (Fig. 4 (a)). Unfor- 
tunately, we were only able to obtain this solution away 
from some small neighbourhood of V = 0. (For h = 0.7, 
the smallest value of the velocity for which we were still 
able to find the solution in question, was V = 0.000283.) 
Whether this branch can be continued to V = Q, remains 
an open question. 



E. Other branches; h < 0.28 

As we have mentioned, for h < 0.28 the branch tp + 
extends all way to V = c where it merges with the zero 
solution. No other solutions can be obtained from the ip+ 
soliton. However, in this case we can obtain new branches 
by continuing the (quiescent) twist soliton, Eq.(p2|) with 

* = c> 

It is convenient to start our description with the mo- 
tionless twist solution with the negative momentum. As 
we move in the direction of positive V, the twist gradu- 
ally transforms into a bound state of two ^--solitons (see 
Fig. 6 (a)). At some V — V max the branch turns back, 
shortly after which, at the point V = V c , the momen- 
tum reaches its maximum and starts decreasing. Adja- 
cent to the turning point is a small range of velocities 
V c < V < Vmax where we have two stable solutions cor- 
responding to each V. If we continue the branch with 
positive-momentum twist solution, also in the direction 
of positive V, the solution gradually transforms into a 



complex of two solitons. The momentum reaches its 
maximum, starts decreasing, then the branch turns back 
in V and we find ourselves approaching the origin on the 
(V,P)-plane (Fig. 6 (a)). As we move towards the origin 
along the 4>(++) or along the V^--) branch, the separa- 
tion between two solitons constituing the corresponding 
complex grows while the imaginary part of the solution 
tends to zero. Similarly to what we had for larger h (sec- 
we conjecture that the separation becomes 
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Fig. 6 (a): The full bifurcation diagram for h — 0.05. 
Thick respectively thin lines depict stable respectively un- 
stable branches. Note a region of stability of the complex 
Additional branches can be generated by employing 
the reflection symmetry V — * —V, P — > —P. (b): The 
corresponding E{P) diagram. 

Fig. 6 (b) shows the corresponding E(P) dependence. 
As in Fig. 4 (b), the cusps mark points of the zero cross- 
ing by the stability eigenvalues. Note that as in the case 
of large h (Fig. 4 (b)), there are solitons with the same 
value of the momentum but different energies. Similarly 
to Fig. 4 (b), the stable branch of if>- has the lowest en- 
ergy; bound states on the stable ?/>t — ► "0(++) branch 
also have lower energies than their counterparts with the 
same P and smaller \V\. However, in the case of the ip+ 
solitons we have an interesting reverse of fortunes: out 
of the two branches with the same P, the stable branch 
is the one with higher energy! 
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V. NONLINEAR STAGE OF INSTABILITY 

In this section we present results of our numerical sim- 
ulations of the full time-dependent nonlinear Schrodinger 
equation d) (with 7 = 0.) The objective was to 
study the nonlinear stage of the development of insta- 
bilities reported in the previous section and to iden- 
tify the attractors emerging as t — > oo. We utilised 
a split-step pseudospectral method, with 2 11 = 2048 
modes on the intervals —40 < X < 40 and —80 < 
X < 80, and with 2 12 = 4096 modes on the inter- 
val (—60,60). The method imposes periodic bound- 
ary conditions ip{L/2,t) = ip(-L/2,t), ip x {L/2,t) = 
1>x(-L/2,t). 

We have simulated the evolution of moving solitons un- 
stable against an oscillatory mode and those with a pos- 
itive, nonoscillatory, eigenvalue in their linearized spec- 
trum. One of our conclusions here is that both types of 
instability give rise to the same asymptotic attractors. 
(This is in agreement with earlier simulations of motion- 
less solitons 

A. The decaying breather 

Depending on the value of the driving strength, the 
initial conditions and the choice of the parameters of the 
numerical scheme, we observed one of the two scenarios. 
In the first scenario the soliton transforms into a bell- 
shaped structure, with a small amplitude and large spa- 
tial width, oscillating approximately as ip ~ e lujt , with 
negative to. This localised solution was previously en- 
countered in numerical simulations of Ref. [fl2f where 
it was termed breather. The amplitude of the breather 
slowly decays with time and the width slowly grows. 

We have detected this scenario for the driving strength 
h = 0.1, with the initial condition in the form of the ip+ 
soliton travelling with the velocity V = 0.05 and with 
V = 0.8. (For both values of the velocity the ip+ soliton 
is unstable against an oscillatory mode.) Unlike earlier 
simulations [|12| which started with the initial condition 
in the form of a quiescent unstable soliton and gave rise 
to a quiescent breather, the breather emerging from a 
travelling soliton has a nonzero speed. 

One may naturally wonder whether the speed of the 
breather will decay to zero or approach a nonzero con- 
stant value as t increases. Our simulations seem to sup- 
port the latter hypothesis. In one run, the speed of the 
breather evolving out of the soliton travelling with the 
initial velocity of V = 0.05, was seen to slowly grow and 
gradually approach the constant value of 0.1. This sim- 
ulation was repeated, with the same parameters of the 
numerical scheme and an initial condition which was only 
different from the previous one due to interpolation er- 
rors of order 10 -6 . In this run the breather was first seen 
to slow down, stop but then start moving in the opposite 
direction with the velocity close to -0.2; see Fig. 7 (a). 



(This remarkable sensitivity to the initial data deserves 
a separate comment; see below.) The velocity of the 
breather evolving out of the V = 0.8 soliton, was tending 
to approximately 2.1. However, for large t the unam- 
bigous interpretation of the numerical data is hindered 
by the growth of the amplitude of the radiation back- 
ground. The radiation waves emitted by the oscillating 
breather re-enter the interval via the periodic boundary 
conditions and at a certain stage their amplitudes become 
comparable with the amplitude of the breather. Conse- 
quently, the constant-velocity motion of the breather may 
have been induced by the interaction with the backround 
radiations. 




(b) _ 40 



Fig. 7 The two types of asymptotic attractors resulting 
from the decay of the unstable steadily travelling solitons: 
(a) the decaying and (b) the growing breather, (a) corre- 
sponds to h = 0.1 and the initial condition in the form of 
the ip + soliton with V = 0.05. In (b), h = 0.05 and the ini- 
tial condition was chosen as the ip_ soliton with V = 0.05. 
In both plots the emerging breather changes, spontaneously, 
its direction of motion. (Note that this happens not as a 
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result of the reflection from the boundary, as the periodic 
boundary conditions are imposed.) 



B. The growing breather 

The decaying breather was detected in simulations on 
the interval (—40,40) with N — 2 11 modes. Changing 
the numerical parameters produced an entirely different 
scenario, however, — for the same value of the control 
parameter in Eq.(|l|) (h — 0.1), and for the same initial 
conditions (V = 0.05 and V = 0.8). 

Namely, we increased the number of the Fourier modes 
to N = 2 12 and the length of the interval first to L = 120 
and then to 160. As in the case of L = 80 and N = 2 , 
in simulations with the new values of iV and L the un- 
stable travelling soliton ip + was seen to transform into a 
bell-shaped structure, oscillating roughly as tp ~ e tu}t . 
However, this time the emerging breather has a posi- 
tive frequency w; its amplitude is large and continues 
to slowly grow, while the width is narrow and keeps on 
decreasing (Fig. 7 (b)). 



This attractor was also observed previously in |12|. It 
was found there that the decaying and growing breather 
coexist. Whether the evolution of the same unstable soli- 
ton settles to one or the other asymptotic attractor, was 
found to depend on the choice of the phase of a small per- 
turbation applied to the initial condition. In our present 
simulations, the perturbation is introduced simply by 
changing the parameters of the numerical scheme. 

We also examined initial conditions in the form 
of translationally unstable solitons, including the ip+~ 
soliton with V = 1.4 for the driving strength h = 0.1 
and the -0_-soliton with initial velocities V = 0.05 and 
V = 1.4, for the driving strength h = 0.05. For each of 
the above three situations the simulations were repeated 
with 2 11 modes on the interval —40 < X < 40, and with 
2 12 modes on the intervals (-60,60) and (-80,80). In 
all nine runs the unstable soliton was seen to evolve into 
the growing breather. (It is quite likely that some other 
choices of the numerical parameters may give rise to the 
decaying breather instead.) 

The velocity of the growing breather may vary during 
its evolution. It can even wander erratically, changing 
the direction of its motion several times, but eventually, 
for t ~ 10 4 or even earlier, the speed of the breather locks 
on to some constant value. Since the amplitudes of ra- 
diation waves are comparable with the amplitude of the 
breather at that stage, this effect can be induced by the 
breather-radiation interactions. 



VI. CONCLUDING REMARKS AND OPEN 
PROBLEMS 

The main result of this paper is the demonstration of 
the existence of wide classes of travelling soliton solu- 



tions of the (undamped) parametrically driven nonlinear 
Schrodinger equation. We established the necessary con- 
ditions under which motionless solitons can be continued 
to nonzero velocities, and, in cases where these condi- 
tions were met, were indeed able to carry out the numer- 
ical continuation. The stability of all resulting branches 
of solutions was examined; oscillatory and translational 
instabilities identified, and the single-soliton stability 
chart compiled on the (h, T^)-plane. The onset values of 
the translational instabilities, obtained numerically, were 
shown to verify the relation dPj dV — predicted by our 
theoretical analysis. As opposed to the case of the soli- 
ton ip- , which undergoes similar transformations for any 
h, the result of the continuation of the has turned 
out to be sensitive to the value of the driving strength. 
We have identified two different transformation scenar- 
ios, one occurring for small and the other one for larger 
h, and described an interesting cross-over from one to the 
other. 

In our analysis we paid a special attention to the tra- 
jectories of linearised eigenvalues on the complex plane. 
Apart from the information on the stability of different 
branches of solutions, the behaviour of the eigenvalues 
can give insight into the supercritical dynamics of soli- 
tons, i.e. dynamics beyond the instability threshold Jl2| . 
The motion of the eigenvalues in the undamped case that 
we are currently concerned with, allows one to even pre- 
dict the asymptotic attractors arising when there is a 
small but nonzero damping in the system jl^]. Rele- 
gating the corresponding bifurcation analysis to future 
publications, here we have restricted ourselves to a se- 
ries of numerical simulations of the time-dependent NLS 
equation (0) (with 7 = 0.) 

It is worth listing, separately, stable solutions obtained 
in this study. First, the quiescent soliton which is 
always unstable with respect to a nonoscillatory mode, 
stabilizes when travels faster than a certain critical veloc- 
ity. The stability boundary satisfies dP/dV — 0. Second, 
the quiscent soliton ip + is stable for h < 0.064 and loses 
its stability to an oscillating soliton for h > 0.064. For 
nonzero V the stability region is shown in Fig. 3. The 
lower boundary of this region corresponds to the onset of 
the oscillatory instability while along the upper boundary 
the soliton becomes unstable with respect to a nonoscil- 
latory mode. The corresponding critical velocity satisfies 
dP/dV = 0. Third, the bound state also displays 

a region of stability for small h — see Fig. 6. 

We conclude this section by pointing out to several 
open questions. In the first place, it would be interest- 
ing to continue the ^'(_T-)-branch in Fig. 4. Will this 
multisoliton complex keep on attaching more solitons 
on its flanks? Another "open-ended" branch {4>(tt)) 
approaches the origin vertically down in Fig. 4 (a); it 
would be interesting to continue this branch as well. The 
striking difference between the bifurcation diagrams for 
h > 0.28 (Fig.4) and h < 0.28 (Fig.6) is also worthy of 
a deeper analysis. Is the bifurcation diagram in Fig.6 
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complete, or there are other multisoliton branches simi- 
lar to those arising for large h (Fig. 4)? The next open 
question concerns the decaying and growing breather so- 
lutions arising as a result of the growth of the instability 
of steadily travelling solitons. If the radiations were pre- 
vented from re-entering the interval of simulation, would 
the velocity of the breather decay to zero or approach 
some nonzero value determined by the initial velocity of 
the unstable soliton? Finally, a challenging problem is 
to find out whether travelling solitons can exist in the 
presence of damping. The effect of dissipation will be, 
of course, to attenuate the soliton. It is not obvious 
whether the spatially uniform parametric pumping is ca- 
pable of compensating this type of losses and sustaining 
the soliton's steady motion. It is fitting to note here that 
the parametrically driven, damped nonlinear Schrodinger 
equation has wide classes of stationary solitons H] , some 
of which may be continuable to nonzero V. 
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